Cluster Monte Carlo Algorithm for the Quantum Rotor Model 
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We propose a highly efficient "worm" like cluster Monte Carlo algorithm for the quantum rotor 
model in the link-current representation. We explicitly prove detailed balance for the new algorithm 
even in the presence of disorder. For the pure quantum rotor model with /i = the new algorithm 
yields high precision estimates for the critical point K c = 0.33305(5) and the correlation length 
exponent v = 0.670(3). For the disordered case, fj, = i ± i, we find v = 1.15(10). 
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What types of insulating, conducting, superconducting 
and more exotic phases occur in two-dimensional systems 
at T = is a topic of considerable current interest. A 
significant amount of theoretical [1-4] and experimen- 
tal [5,6] effort has focused on bosonic systems where a 
superconductor to insulator transition is known to occur. 
In agreement with most experiments [5], it was under 
quite general conditions shown [1,2] that a transition can 
occur directly between the superconducting and insulat- 
ing states. However, more recently, it has been suggested 
that an exotic metallic phase also is possible [4,6]. In this 
context precise numerical results would be very valuable 
and in the present paper we propose a new, very efficient 
cluster Monte Carlo algorithm for this purpose, allowing 
us to significantly improve previous results. In particular 
we show that the inequality [7] v > 2 j d is not violated in 
the presence of disorder, resolving contradictions in pre- 
vious work. The high precision of the algorithm should 
allow for precise calculations of transport properties of 
quantum rotor models studied theoretically in [1,4]. The 
ideas presented here could be useful for the study of clas- 
sical spin systems [19]. 

Low-dimensional bosonic systems are often described 
in terms of the (disordered) boson Hubbard model: 
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U is the on-site repulsion, to the hopping strength, /i r the 
chemical potential varying uniformly in space between 
/.t ± A and h r = $J$ r is the number operator. If we set 
$ r = | < I , r |e* 9r and integrate out amplitude fluctuations, 
i?bH becomes equivalent to the quantum rotor model [8] : 
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Here, r the phase of the quantum rotor, t the renor- 



malized hopping strength and i- 



The quantum 



rotor model describes a wide range of phase transitions 
dominated by phase-fluctuations and it is well known [8] 
that an equivalent classical model exists where the Hamil- 
tonian is written in terms of currents defined on the links 



of a lattice, J — (J x , J v , J T ). These link-current vari- 
ables describe the "relativistic" bosonic current which 
should be divergenceless, V • J = 0. In terms of these 
variables the classical (2+l)D Hamiltonian can be writ- 
ten as follows [81: 
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E' denotes a summation over configurations with V • J = 
0. Varying K corresponds to changing the ratio t/U in 
the quantum rotor model. The quantum rotor model has 
been extensively studied [8-10] in this representation, but 
a number of conclusions can be questioned due to severe 
finite-size effects. For notational convenience it is useful 
to slightly enlarge the definition of the link-currents at 
a given site in the following way: At each site (r, r) in 
the lattice we define six surrounding link variables JZ, T \ 
where a runs over ±x, ±y, ±r. Note that, with this no- 
tation JZ, x y T j and J( x _i y r ) is the same variable, with 
equivalent relations in the other directions. The diver- 
genceless constraint at the site (r, t) can then be written: 

J-x + j-y + j-t _ jx _j_ jy + jt^ g0 that thc gum Q £ the 

incoming and outgoing currents are equal. Conventional 
Monte Carlo updates [8] on this model consists of updat- 
ing simultaneously four link variables. Global moves, up- 
dating a whole line of link variables thus allowing particle 
and winding numbers to fluctuate, are added to ensure 
ergodicity, but the acceptance ratio for these moves be- 
comes exponentially small for large lattice sizes. Here we 
will describe a way to construct a worm-like algorithm to 
perform non local moves for this model. 

The cluster algorithm [11-13] we propose is similar in 
spirit to worm algorithms [14,15] in the sense that we up- 
date the link-currents by moving a "worm" through the 
lattice visiting the sites s% = (r^r^). The links through 
which thc worm pass are updated during its construction. 
At a given site, the links with er equal to x,y,T are called 
outgoing links and those with a equal to — x, —y, — r in- 
coming links. When the worm is moving through the lat- 
tice the currents Jf. are updated in the following manner: 
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if the worm is leaving the site s, along an outgoing link 
we increment the corresponding current: 



K J'Z = Jl+h tr = x,y,T. 



(3) 



If the worm is leaving the site Sj along an incoming link 
we decrement the corresponding current: 



j£ = •/£-!, a = -s, -y, -t. 



(4) 



The construction of the worm starts with the choice of a 
random initial site s\ — (ri, ti) in the space-time lattice. 
Then the algorithm can be decomposed in two steps, (i) 
The worm moves to one of the 6 neighboring sites. To 
decide which direction to go from a site Si — (rj,Tj), we 
calculate for all directions a = ±x, ±y, ±r weights, A* , 
according to local detailed balance. A good choice is: 

A° t = min(l,exp(-A££/JO). A££ = E'° - ££. (5) 

Here E a s . = |(Jf.) 2 — /Un^^CT,±r is the contribution to 
the total energy from the link J". , before the worm moves 
through it. E'° is the energy contribution with J% re- 
placed by J 1 ". By normalizing the A°. 's we define the 
probabilities:* p°. = A a s jN s%1 where N s% = ^2 a A" Si . A 
direction a is then chosen according to these probabili- 
ties, (ii) Once a is chosen, we update the corresponding 
J£ according to the above rules, Eq. (3) and (4), and 
extend the worm to the new lattice site Sj+i. (i) and 
(ii) are then repeated until the worm passes through the 
initial site where Sj+i = s\. Finally, in order to satisfy 
detailed balance we have to erase the worm with a prob- 
ability determined in the following way. If iV(worm) and 
iV(no worm) are the normalization of the probabilities 
at the initial site s\ with and without the worm present, 
then we erase the constructed worm with a probability 
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Under most conditions this probability is very small. Sev- 
eral points are noteworthy about this algorithm. First of 
all, the configurations generated during the construction 
of the worm are not valid (V • J ^ 0). However, once 
the construction of the worm is finished and the path of 
the worm closed, the divergenceless constraint is satis- 
fied. Secondly, when the worm moves through the lattice 
it may pass many times through the same link and cross 
itself before it reaches the initial site where the construc- 
tion terminates. Hence, it is crucial that the current vari- 
ables are updated during the construction of the worm. 
Finally, at each step i in the construction of the worm it 
is likely that the worm at the site Sj will partially "erase" 
itself by choosing to go back to the site Sj_i visited im- 
mediately before, thereby "bouncing" off the site Sj. 

Now we turn to the proof of detailed balance for the 
algorithm. Let us consider the case where the worm, w, 
visits the sites {si ... sat} where s\ is the initial site. The 



worm then goes through the corresponding link variables 
{h . . -In}, with U connecting S{ and Sj+i. Note that 
sjv is the last site visited before the worm reaches s\. 
Hence, sn and s\ are connected by the link In- The total 
probability for constructing the worm w is then given by: 

P w = P S1 (1-P*) Ilili A sJ N s> ■ The indcx a denotes the 
direction needed to go from Sj to Sj+i , P Sl is the probabil- 
ity for choosing site s\ as the starting point and is the 
probability for erasing the worm after construction. If the 
worm w has been accepted we have to consider the prob- 
ability for reversing the move. That is, we consider the 
probability for constructing an anti-worm w annihilating 
the worm to. We have: P^ = P Sl (1 - P|) n^i A %JN- Si . 
Here, the indcx a denotes the direction needed to go from 
Si to Sj+i. Note that, in this case the sites are visited in 
the opposite order, s\ = S\,S2 — Sn,---,Sn — S2, in 
general Sj = SN- i+2 (i ^ 1). Note also that, Sj and 
Sj+i are connected by the link k = In-i+i with sn and 
s~i connected by In = h- With this notation we see 
that, Si and sjv-i+i = Sj+i are connected by the link 
variable ij. Let us now consider the case where both 
of the worms w and w have reached the site Si differ- 
ent from the starting site s\. Since we are updating the 
link variables during the construction of the worm and 
since we are always considering moving the worm in all 
six directions, we have N Si = N SN _ i+2=Si (i ^ 1). Fur- 
thermore, AZ and A% _„ only depend on the link 
variable U connecting the sites si and sn-i+i and we 
see that: A"JA% S _^ = cxp(-AE°jK), i = 1 . . . N. 
Hence, since P Sl = P Sl , we find: 
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where Ai?Tot is the total energy difference between a con- 
figuration with and without the worm w present. Now 
we consider P e = 1 — min(l, iV Sl (no worm)/iV Sl (worm)). 
Here, iV Sl = N Sl (no worm) is equal to N Sl (anti-worm) and 
iV Sl (no anti-worm)=7V Sl is equal to A f Sl (worm). Hence, 
we find for the probability to erase the worm P^ = 
1 - min(l, N S1 /N S1 ) and P% = 1 - min(l, N s jN Sl ) for 
erasing the anti-worm. With this choice of P e we satisfy 
detailed balance since: -p^- = exp(— A.E^ ot /K). Ergodic- 
ity is simply proven as the worm can perform local loops 
and wind around the lattice in any direction, as in the 
conventional algorithm. 
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FIG. 1. Autocorrelation times versus lattice size for the 
conventional and worm algorithm for p = at K = 0.333. 
The dashed lines indicate power-law fits and the solid line an 
exponential fit in L. 
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FIG. 2. Lp versus K for different lattice sizes, for p = 0. 
All curves cross at the critical point K c = 0.33305(5) with 
Lp\x=K c = 0.495(5). Inset: Ldp/dK at K c versus L. The 
dashed line indicates a fit yielding an exponent v — 0.670(3). 



To demonstrate the efficiency of the proposed algo- 
rithm we have calculated auto-correlation times for dif- 
ferent lattice sizes for the worm algorithm and the con- 
ventional algorithm. For an observable O we define the 
auto-correlation function and the auto-correlation time 
to in the usual manner [16]: 

Here, t is the Monte Carlo time measured in Monte 
Carlo sweeps (MCS), with 1 MCS corresponding to L d 
attempted updates. The auto-correlation function is cal- 
culated from simulations with 10 8 MCS, and to obtain 
the best estimate of to we always fit to the indicated 
double-exponential form with n C tq. To make a fair 
comparison of to for the two algorithms, one custom- 
arily [12,16] multiplies To for the worm algorithm by 
N / (I), with (l) the mean number of links in a worm and 
N = 3L 3 . With this rescaling we show in Fig. 1 the auto- 
correlation times, t p for the stiffness (see exact definition 
below) at p = for both algorithms. The calculations 
have been performed on cubic lattices at K = 0.333, very 
near previous estimates of the critical point [9]. For the 
worm algorithm we also show the auto-correlation time 
for the energy, te, which is almost identical to t p . The 
auto-correlation times increase dramatically with system 
size for the conventional algorithm where as they remain 
very small (of the order of 2-3 MCS per link) for the worm 
algorithm. If we fit the L dependence of t p ~ L ZMC 
with a power law, we obtain an auto-correlation expo- 
nent Zmc larger than 4 for the conventional algorithm. 
For the conventional algorithm it is likely that r p is di- 
verging exponentially with L since p is solely determined 
by global updates for which the acceptance probability 
decreases exponentially with L. For the worm algorithm 
we find a very small zmc ~ 0.3. 



We now present results for the model Eq. (2) at p, = 0. 
There, the model is expected to undergo a transition in 
the (2+l)D XY universality class [1,9] from a superfluid 
into a Mott insulating phase with a dynamical critical ex- 
ponent z = 1. The different phases can be distinguished 
by calculating the stiffness defined as [8] : 

P=Tj2^ J (r,r)) 2 )- (9) 

Since we expect z = 1, we use L T , the system size in the 
third direction, equal to L. To obtain the K dependence 
of p we have used reweighting techniques [17] on large 
runs (of the order of 10 8 MCS) at K = 0.333. The er- 
ror bars are determined using jackknife techniques [16]. 
Using finite-size scaling relations, the quantity pL z is ex- 
pected to be independent of system size at the critical 
point [8], K c . Moreover, L z dp/dK is expected to di- 
verge at K c as L x l v where v is the correlation length 
exponent. We have explicitly calculated this quantity by 
evaluating the thermodynamic derivative of p with re- 
spect to the coupling K: dp/dK = ((pE) - {p){E))/K 2 . 
In Fig. 2, we show Lp versus K for different lattice 
sizes. From the crossing of the curves we can deter- 
mine K c = 0.33305(5) to a much higher precision than 
was possible using the conventional algorithm on much 
smaller systems [8-10]. Since all the curves cross in a 
single point our results are completely consistent with a 
dynamical exponent z = 1, as expected [1]. In the in- 
set of Fig. 2 is shown the size dependence of L z dp/dK 
at K c on a log-log scale. We fit this curve to a power- 
law AL 1 / 1 ' and obtain v = 0.670(3), in perfect agreement 
with estimates for the three-dimensional XY universal- 
ity class [18]. Preliminary results [19] for the generic 
transition at pi = \ show pronounced finite-size effects 
questioning previous work [10]. 
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FIG. 3. (A) [L 2 p] av versus K for different lattice sizes, for 
fj, = i ± |. All curves cross at the critical point K c — 0.246(1) 
with [L 2 p] !IV \k=k c =0.12(1). Inset: [L 2 dp/dK] av versus L for 
different K . The solid line indicates a power-law fit yielding 
an exponent v = 1.15(10). (B) Scaling plot of L 2 p(L,L T ) at 
K c = 0.246. 



We also simulated the model Eq. (2) with disorder 
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In this case the transition is between 



for fi ■- 2 _ 2 . 
a supcrfluid and an insulating boseglass phase. Scal- 
ing theory [1] predicts a second order transition with 
dynamical exponent z = 2. Hence, we use lattices of 
size L x L x aL 2 where a = L T /L 2 is the aspect ratio. 
Previous work [8], limited to L < 10, have determined 
K c = 0.248 ± 0.002. Estimates for the correlation length 
exponent [8,10] yielded v — 0.9 ±0.1 almost violating the 
inequality [7] v > 2/d. From the results shown in Fig. 3 
(A), obtained with the cluster algorithm, it is clear that 
K c in fact is at a slightly lower value K c — 0.246(1), al- 
though the crossing of L = 6, 8 occurs at K = 0.248. The 
disorder average, [-] av , has been performed over 50,000 
samples with 10 5 MCS per sample. The more precise 
value for K c significantly changes estimates of v. The 
inset in Fig. 3 (A) shows [L 2 dp/dK] av versus L, which 
at K c yields v — 1.15(10) now largely satisfying the in- 
equality v > 2/d. The results in Fig 3 (A) are clearly 
consistent with z = 2. In Fig. 3 (B) we show results for 
L 2 p(L,L T ) versus L T /L 2 at K c . Standard scaling the- 
ory [20] predicts that this should be a universal function 
of a if z — 2. Our results nicely confirm this. The values 
of exponents are in good agreement with the analytical 
estimates in ref. [21]. 



In conclusion, we have introduced a worm algorithm 
for the quantum rotor model. For the link-current rep- 
resentation of the quantum rotor model the proposed 
algorithm is exponentially more efficient than conven- 
tional algorithms and performs at par with the Wolff 
algorithm [12] for the classical 3D XY model. Most note- 
worthy, the algorithm performs exceptionally well on dis- 
ordered systems. We have also successfully adapted it to 
the study of systems with longer range interactions as 
well as classical Ising models [19]. 
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